#ifndef _POLYNOMIAL_H_
#define _POLYNOMIAL_H_

#include "Vec.h"
#include "Numtool.h"
#include <vector>

#define Max(a,b) ((a)>(b))?(a):(b)

template<int Order, class CoefType> class Polynomial;

template<int Order, class CoefType>
std::ostream &operator<<(std::ostream & os ,const Polynomial<Order,CoefType> &p);
template<int Order1, int Order2, class CT> Polynomial<Max(Order1,Order2),CT>
operator+(const Polynomial<Order1,CT>& p1,const Polynomial<Order2,CT>& p2);
template<int Order1, int Order2, class CT> Polynomial<Max(Order1,Order2),CT>
operator-(const Polynomial<Order1,CT>& p1,const Polynomial<Order2,CT>& p2);
template<int Order1, int Order2, class CT> Polynomial<Order1+Order2-1,CT>
operator*(const Polynomial<Order1,CT>& p1,const Polynomial<Order2,CT>& p2);
template<int Order, class CoefType> Polynomial<Order,CoefType>
operator*(const Polynomial<Order,CoefType>& p, const double x);
template<int Order, class CoefType> Polynomial<Order,CoefType>
operator/(const Polynomial<Order,CoefType>& p, const double x);



template<int Order, class CoefType>
class Polynomial{
protected:
  CoefType coe_[Order];
public:
  Polynomial(const CoefType& coe = CoefType());
  Polynomial(const std::vector<CoefType>& vcoe);
  Polynomial(std::initializer_list<CoefType> l);
  ~Polynomial() = default;
  const CoefType& operator[] (const int k) const;
  CoefType& operator[] (const int k);
  template<int Order1, int Order2, class CT> friend
  Polynomial<Order1,CT> adjust(const Polynomial<Order2,CT>& p);
  template<class T> CoefType operator() (const T& x) const;
  Polynomial<Order+1,CoefType> integral() const;
  Polynomial<Max(1,Order-1),CoefType> der() const;
  template<class T> Polynomial<Order,CoefType> translate(const T& x0) const;
  friend std::ostream &operator<< <>(std::ostream &os, const
  Polynomial<Order,CoefType> &p);
  template<int Order1, int Order2, class CT> friend Polynomial<Max(Order1,Order2),CT> operator+(const
  Polynomial<Order1,CT>& p1, const Polynomial<Order2,CT>& p2);
  template<int Order1, int Order2, class CT> friend Polynomial<Max(Order1,Order2),CT> operator-(const
  Polynomial<Order1,CT>& p1, const Polynomial<Order2,CT>& p2);
  template<int Order1, int Order2, class CT> friend Polynomial<Order1+Order2-1,CT> operator*(const
  Polynomial<Order1,CT>& p1, const Polynomial<Order2,CT>& p2);
  friend Polynomial<Order,CoefType> operator* <> (const
  Polynomial<Order,CoefType>& p, const double x);
  friend Polynomial<Order,CoefType> operator/ <> (const
  Polynomial<Order,CoefType>& p, const double x);
  
};


template<int Order, class CoefType>
inline Polynomial<Order,CoefType>::Polynomial(const CoefType& coe){
  coe_[Order-1] = coe;
  for (int i = 0 ; i < Order-1 ; i++)
    coe_[i] = CoefType();
}

template<int Order, class CoefType>
inline Polynomial<Order,CoefType>::Polynomial(const std::vector<CoefType>& vcoe){
  const int N = vcoe.size();
  assert(Order == N);
  for (int i = 0 ; i < Order ; i++)
    coe_[i] = vcoe[i];
}

template<int Order, class CoefType>
inline Polynomial<Order,CoefType>::Polynomial(std::initializer_list<CoefType> l){
  auto j = l.begin();
  for(int i = 0; i < Order ; i++)
      coe_[i] = *j++;
}

template<int Order, class CoefType>
inline const CoefType& Polynomial<Order,CoefType>::operator[] (const int k) const{
  return coe_[k];
}

template<int Order, class CoefType>
inline CoefType& Polynomial<Order,CoefType>::operator[] (const int k){
  return coe_[k];
}


template<int Order, class CoefType> template<class T>
inline CoefType Polynomial<Order,CoefType>::operator() (const T& x) const{
  CoefType res = coe_[Order-1];
  T basis = static_cast<T>(1.0);
  for (int i = Order-2; i >= 0 ; i--){
    basis = basis*x;
    res = res + coe_[i]*basis;
  }
  return res;
}


template<int Order1, int Order2, class CT>
Polynomial<Order1,CT> adjust(const Polynomial<Order2,CT>& p){
  Polynomial<Order1,CT> res;
  if (Order1 >= Order2){
    for (int i = 1 ; i <= Order2 ; i++)
      res[Order1-i] = p[Order2-i];
  }
  else{
    for (int i = 1 ; i <= Order1 ; i++)
      res[Order1-i] = p[Order2-i];
  }
  return res;
}

template<int Order, class CoefType>
inline Polynomial<Order+1,CoefType> Polynomial<Order,CoefType>::integral() const{
  Polynomial<Order+1,CoefType> res;
  for (int i = Order - 1 ; i >= 0 ; i--)
    res[i] = coe_[i]/(Order - i);
  return res;
}


template<int Order, class CoefType>
inline Polynomial<Max(1,Order-1),CoefType> Polynomial<Order,CoefType>::der()
const{ 
   Polynomial<Max(1,Order-1),CoefType> res;
   for (int i = Order - 2 ; i >= 0 ; i--)
      res[i] = coe_[i]*(Order - i - 1);
   return res;
}

template<int Order, class CoefType> template<class T>
Polynomial<Order,CoefType> Polynomial<Order,CoefType>::translate(const T& x0) const{
  Polynomial<Order,CoefType> res;
  Polynomial<Order,CoefType> tmp = *this;
  res[Order - 1] = (*this)(x0);
  for (int i = Order - 2 ; i >= 0 ; i--){
    tmp = adjust<Order,Order-1>(tmp.der());
    res[i] = tmp(x0)/factorial(Order - i - 1);
  }
  return res;
}


template<int Order, class CoefType>
inline std::ostream& operator<<(std::ostream &os, const
Polynomial<Order,CoefType> &p){
  const char c = 't';
  for (int i = Order; i > 2 ; i--){
    os << p.coe_[Order-i] << c << "^" << i-1 << "+";
  }
  if (Order > 1)
    os << p.coe_[Order-2] << c << "+";
  if (Order > 0)
    os << p.coe_[Order-1];
  else
    os << CoefType();
  return os;
}


template<int Order1, int Order2, class CT>
inline Polynomial<Max(Order1,Order2),CT> operator+(const
Polynomial<Order1,CT>& p1, const Polynomial<Order2,CT>& p2){
  if (Order1 >= Order2){ 
    Polynomial<Max(Order1,Order2),CT> res;
    const int delta = Order1 - Order2;
    int i;
    for (i = Order1-1 ; i >= delta ; i--)
      res[i] = p1[i] + p2[i-delta];
    for(; i >= 0 ; i--)
      res[i] = p1[i];
    return res;    
  }
  else
    return p2+p1;    
}


template<int Order1, int Order2, class CT>
inline Polynomial<Max(Order1,Order2),CT> operator-(const
Polynomial<Order1,CT>& p1, const Polynomial<Order2,CT>& p2){
  if (Order1 >= Order2){
    Polynomial<Max(Order1,Order2),CT> res;
    const int delta = Order1 - Order2;
    int i;
    for (i = Order1-1 ; i >= delta ; i--)
      res[i] = p1[i] - p2[i-delta];
    for(; i >= 0 ; i--)
      res[i] = p1[i];
    return res;    
  }
  else{
    Polynomial<Max(Order1,Order2),CT> res;
    const int delta = Order2 - Order1;
    int i;
    for (i = Order2-1 ; i >= delta ; i--)
      res[i] = p1[i-delta] - p2[i];
    for(; i >= 0 ; i--)
      res[i] = -p2[i];
    return res;     
    }
}


template<int Order1, int Order2, class CT>
inline Polynomial<Order1+Order2-1,CT> operator*(const
Polynomial<Order1,CT>& p1, const Polynomial<Order2,CT>& p2){
  Polynomial<Order1+Order2-1,CT> res;
  for (int i = 0 ; i < Order1+Order2-1 ; i++)
    for (int j =std::max(0,i-Order2+1); j <= std::min(i,Order1-1) ; j++)
      res[i] = res[i] + p1.coe_[j] * p2.coe_[i-j];
  return res;
}


template<int Order, class CoefType>
inline Polynomial<Order,CoefType> operator*(const
Polynomial<Order,CoefType>& p, const double x){
  Polynomial<Order,CoefType> res;
  for (int i = 0 ; i < Order ; i++)
    res[i] = p.coe_[i] * x;
  return res;
}

template<int Order, class CoefType>
inline Polynomial<Order,CoefType> operator/(const
Polynomial<Order,CoefType>& p, const double x){
  Polynomial<Order,CoefType> res;
  for (int i = 0 ; i < Order ; i++)
    res[i] = p.coe_[i] / x;
  return res;
}



#endif // _POLYNOMIAL_H_
